#Nature's Kidneys: the Role of Wetland Reserve Easements in Restoring Water Quality
#Nicole Karwowski and Marin Skidmore
#Replication File
#2025

#Extract PRISM data to each huc12

#load libraries
library(ggplot2)
library(raster)
library(sf)
library(exactextractr)
library(prism)
library(tidyverse)
library(mapview)

#set wd
setwd("C:/Users/16308/OneDrive - Montana State University/UW_Madison/Water quality/Replication Package")

################################################################################
#Precipitation

#open prism files
prism_set_dl_dir("C:/Users/16308/OneDrive - Montana State University/UW_Madison/Water quality/Replication Package/Water quality/Data/PRISM/Raw/ppt_stack")

#download files (only need to do this once)
#snapd data goes from '1980-01-01', '2018-12-31'
get_prism_monthlys(type = "ppt", year = 1980:2018, mon = 1:12, keepZip = FALSE)

#create raster stack
stack<-pd_stack(prism_archive_subset(type="ppt", temp_period="monthly", 
                                     year = 1980:2018, mon = 1:12))
crs(stack)

#open watershed files
huc12<-st_read("Data/WatershedBoundaries/Raw/USGS_WBD/USGS_WBD.shp")

#just keep huc12 and geometry
huc12<-subset(huc12, select=c(huc12, geometry))

#turn into same crs as prism data
huc12<-st_transform(huc12, crs(stack))
crs(stack)
crs(huc12)

#extract
#try to extract
ppt_huc12_stack<-exact_extract(stack, huc12, 'mean', progress=T)
names(ppt_huc12_stack)
list(ppt_huc12_stack)

names(ppt_huc12_stack)<-substring(names(ppt_huc12_stack), 29,34)

ppt_huc12_stack$huc12<-huc12$huc12
colnames(ppt_huc12_stack)

summary(ppt_huc12_stack)

#reshape
ppt<-pivot_longer(ppt_huc12_stack, cols=1:468, names_to = "YearMonth", values_to = "Ppt")

summary(ppt)

#lets save it
save(ppt, file="Data/PRISM/Processed/ppt_huc12_monthly_1980_2018.RData")


################################################################################
#Mean temperature 

#set wd
prism_set_dl_dir("C:/Users/16308/OneDrive - Montana State University/UW_Madison/Water quality/Replication Package/Water quality/Data/PRISM/Raw/temp_stack")

#download files (only need to do this once)
#snapd data goes from '1980-01-01', '2018-12-31'

get_prism_monthlys(type = "tmean", year = 1980:2018, mon = 1:12, keepZip = FALSE)

#create raster stack
stack<-pd_stack(prism_archive_subset(type="tmean", temp_period="monthly", 
                                     year = 1980:2018, mon = 1:12))
crs(stack)

#open watershed files
huc12<-st_read("Data/WatershedBoundaries/Raw/USGS_WBD/USGS_WBD.shp")

#just keep huc12 and geometry
huc12<-subset(huc12, select=c(huc12, geometry))

#turn into same crs as prism data
huc12<-st_transform(huc12, crs(stack))
crs(stack)
crs(huc12)

#extract
tmp_huc12_stack<-exact_extract(stack, huc12, 'mean', progress=T)
names(tmp_huc12_stack)
list(tmp_huc12_stack)

names(tmp_huc12_stack)<-substring(names(tmp_huc12_stack), 31,36)

tmp_huc12_stack$huc12<-huc12$huc12
colnames(tmp_huc12_stack)

summary(tmp_huc12_stack)

#reshape 
tmp<-pivot_longer(tmp_huc12_stack, cols=1:468, names_to = "YearMonth", values_to = "Temp")

summary(tmp)
table(tmp$YearMonth)

#lets save it
save(tmp, file="Data/PRISM/Processed/tmp_huc12_monthly_1980_2018.RData")

